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Path integral Monte Carlo (PIMC) simulations of liquid helium at negative pressure have been car- 
ried out for a temperature range from the critical temperature to below the superfluid transition. We 
have calculated the temperature dependence of the spinodal line as well as the pressure dependence 
of the isothermal sound velocity in the region of the spinodal. We discuss the slope of the superfluid 
transition line and the shape of the dispersion curve at negative pressures. 
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I. INTRODUCTION 

Consider a liquid that is quenched from above the 
liquid-gas coexistence line to a point below, which will be 
in either the metastable or unstable region, see Fig. [j]. As 
it approaches equilibrium, the system will phase separate 
to a state of positive pressure with coexisting vapor and 
liquid phases of densities p v and pi respectively. Below 
this coexistence line is the spinodal line, which delineates 
the metastable phase from the unstable phase. The spin- 
odal line is the locus of points where the speed of sound 
vanishes, 7714c 2 = dP/dp where 7714 is the mass of a He 
particle. At the spinodal there is no energy barrier to nu- 
cleation and phase separation. If the temperature is low 
enough the pressure at the liquid spinodal, at a density 
of psi, may be negative. Once the spinodal pressure be- 
comes negative, the lowest pressure the system can attain 
is the spinodal pressure. 

Direct measurement of the spinodal pressure of liq- 
uid helium by homogeneous nucleation is experimentally 
difficult. In a driven system- .such as the pressure oscil- 
lation experiments of Maristm and Balibarffa, negative 
pressures are only achieved for .a finite duration. TIicl 
presence of objects like vorticesB, electronsaQ or bothEl 
lower the nucleation energy barrier. Measurements of the 
cavitation pressure are higher than the spinodal pressure 
because of quantum tunneling^ or thermal activationL2l 
over the barrier, but are consistent with predictions of 
the nucleation energy barrier, attempt frequency and the 
spinodal pressure. In this paper microscopic PIMC sim- 
ulations of liquid helium at negative pressure and finite 
temperature will be presented. The calculated temper- 
ature and density dependence of the .spwipdal line com- 
pare favorably with other calculationsUdlZI. We estimate 
that the superfluid transition can be extended to negative 
pressures. Finally, properties of the excitation spectrum 
are found to be consistent with conjectur.es based on the 
quasi-particle model of superfluid hcliun£3. 



II. METHOD 

At zero temperature, density, functionate 
and microscopic/phenomenologicalEjO calculations and 
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FIG. 1. Isotherms of pressure versus density for a generic 
liquid are shown for temperatures at or less than the critical 
temperature. The spinodal line and the coexistence line form 
the lower and upper bounds of the shaded metastable region. 
Locations of the saturated vapor and liquid densities, p v and 
pi, and the spinodal vapor and liquid densities, p sv and p s i, 
for a given isotherm are indicated. 



quantum Monte CarlcOEj simulations of liquid helium 
at negative pressure have produced comparable values 
of the spinodal pressure and density. At finite tempera- 
ture, a density functionaO developed for the liquid-gas 
interface has been used to successfully examine helium at 
negative pressures. Microscopic finite temperature sim- 
ulations at negative pressures have not been performed 
until now. 

It has been well documented that PIMC can provide 
accurate thermodynamic and some dynamical properties 
of quantum systems in equilibrium .. a.nd has been espe- 
cially successful with liquid hcliurnEallj. The basis for the 
path integral method is the evaluation of the many par- 
ticle density matrix, p — exp(— (3TL). The Hamiltonian is 
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assumed to be 



w = -a£ v ' + 5>^-)' (i) 

where an accurate two-body interaction v is used. 
From the density matrix, expectation values of observ- 
ables can be determined, 

(O) = Z- 1 J dRdR'p(R, R'){R'\0\R), (2) 

with the partition function Z , 



dRp(R, R), 



(3) 



and where the matrix elements of the density matrix in 
position basis is 



p(R,R') = (R\e' fm \R'}. 



(4) 



R is the 3iV state of the system, R = {r\,r%, ..Tjv}- Us- 
ing the product property, the density matrix can be ex- 
pressed as 



-pn _ 



tH\M 



(5) 



where Mr = [3 = 1/ksT and r is the imaginary time 
step. By expressing the density matrix at T as a product 
of density matrices at a higher temperature of MT, Eq. || 
and Eq. |^ can be accurately evaluated by Monte Carlo. 
In our simulation, the canonical ensemble is used, with 
a fixed number of particles N, simulation volume and 
temperature T. For 4 He, A = 6.05961 K A 2 , and during 
the simulations r = 0.0125 K" 1 and N = 64. PIMC in- 
corporates Bose statistics, necessary for the modeling of 
superfluid 4 He, by allowing the permutation of pacticle 
paths. Superfluidity manifests itself as the winding^ of 
the particle paths across the simulation cell when peri- 
odic boundaries are present. 

An advantage of a finite system is its ability to explore 
systems at negative pressures. A disadvantage is that fi- 
nite size effects are present, which can be appreciable for 
a small system and can be difficult to correct. The liquid- 
gas interface limits how large a system can be simulated 
due to the surface energy of a phase separated system. 
Because the simulation is in equilibrium it will phase sep- 
arate if it is thermodynamically favorable. Choosing a 
system size that is on the order of the liquid-gas inter- 
face width prevents the system from formingj-a stable 
coexistence. From previous PIMC simulationslla it was 
found that a 64 particle system is large enough to pro- 
vide accurate bulk properties. To check that this holds 
true near the spinodal, the pressure near the spinodal 
was examined for several systems from 8 to 256 particles 
at a temperature of 2.0 K. In the range of 32 to 128, the 
pressure remained negative and varied little with parti- 
cle size, and S(k) < 1 for small k (see below). A system 
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FIG. 2. The pair correlation function and the static struc- 
ture function are shown for 3 densities at a temperature of 
2.0 K: 0.006 g/cc, 0.073 g/cc and 0.153 g/cc. 



of 64 particles was found to be optimal for minimizing 
finite size effects r^ftd keeping the fluid homogeneous at 
negative pressureEZI. 

To differentiate between liquid and gas phases, the 
static structure function, S(k), and the pair distribu- 
tion function, g(r), were examined at three representative 
densities. Information about fluctuations into the coex- 
isting liquid-gas phase, which we wish to avoid, is present 
in both S(k) and g(r). In Fig. || the two functions are 
shown for three densities and a temperature of 2.0 K. 
The high density liquid phase and the low density gas 
phase show the typical features in S(k) and g(r). The 
intermediate density, which is lower than the liquid spin- 
odal density at that temperature, begins to show signs of 
phase separation but no large scale features (S(k) 3> 1 
for small k) are present. As the density is varied from 
liquid to gas the functions should go smoothly from one 
phase to the other, for a finite system. 



III. ANALYSIS 

A. Spinodal Line 

The PIMC density-pressure isotherms are plotted in 
Fig. [| To determine the location of the liquid spinodal 
line, a cubic polynomial in density was fit to each density- 
pressure isotherm in the region of the liquid spinodal. As 
will be shown below, the cube of the isothermal speed of 
sound is seen to vary linearly with pressure and from this 
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FIG. 3. Spinodal pressure from PIMC simulations are 
compared to other calculations and estimations of the spin- 
odal pressure. Increasing temperature runs from bottom to 
top. Two isotherms from PIMC simulations are shown as 
well, with dotted lines to guide the eye. 



FIG. 4. Temperature dependence of the spinodal pres- 
sures and the superfluid transition at negative pressures are 
shown. At low temperatures the spinodal pressure is insensi- 
tive to temperature while at higher temperatures the behavior 
is linear. The upper solid lines form usual phase diagram. 



is can be shown that P ~ P s ; + a(p — p s i) 3 near the liq- 
uid spinodalEl In Table | the liquid spinodal pressure 
and density for various temperatures are listed and are 
plotted in Fig. ||. These points are indicated by the filled 
circles. The x 2 goodness of fit is on the order of 1 for low 
temperatures and on the order of 2 at the highest tem- 
peratures. Simulations of the gas phase were not done 
for all temperatures and so information about the gas 
spinodal and the liquid-gas coexistence densities is not 

available. 

For comparison— density functionaliil'liS, extrapo- 
lated pcperimentaltj, optimized hypernetted chainEH and 
QMOi-3 values are shown in Fig. |. The extrapolated 
values depend on the assumption that the isothermal 
speed of sound has the following power law dependence, 
c 3 ~ P — P s i , where P s i is the liquid spinodal pressure. 
The density functional values are the result of a param- 



TABLE I. 4 He liquid spinodal values. 



T{K) 


Psi (g/cc) 


P sl (bars) 


0.50 


0.1070(51) 


-8.99(3) 


1.00 


0.1077(45) 


-8.93(3) 


1.54 


0.1097(38) 


-8.62(3) 


1.74 


0.1110(51) 


-8.28(4) 


2.00 


0.1123(45) 


-7.77(3) 


2.50 


0.1136(51) 


-6.27(3) 


3.00 


0.1123(51) 


-4.52(1) 


3.64 


0.110(14) 


-2.68(6) 


4.00 


0.107(28) 


-1.51(1) 



eterized density functional that reproduces the experi- 
mental bulk properties but with the Bose condensation 
added post hoc. While there is agreement with the tem- 
perature dependence of the liquid spinodal pressure, the 
PIMC values of the liquid spinodal density tend to be 
higher than the density functional and extrapolated val- 
ues. Linearly extrapolating the PIMC data to K yields 
a spinodal density very close to the QMC, density func- 
tional and hypernetted chain values. The extrapolated 
K PIMC spinodal pressure is slightly larger than the 
other K results. A possible reason for this discrepancy, 
at least with the QMC value, i^_the use of different ver- 
sions of the Aziz pair potentialEi 

In Fig. H the temperature dependence of the liquid 
spinodal pressure, as calculated by PIMC, is plotted. For 
T in the neighborhood of the critical temperature T c , the 
spinodal pressure behaves similar to that of a Van der 
Waals gas with the spinodal pressure decreasing mono- 
tonically with decreasing temperature. In this region we 
see good agreement of the spinodal pressures. At temper- 
atures near zero, the spinodal pressure is approximately 
independent of temperature, being slightly higher than 
other values when extrapolated to K . This flat region 
may be understood in terms of the quasi-particle picture 
of liquid heliumli 2 ] and will be discussed below. At tem- 
peratures near 2 K there is a distinct transition between 
the two regions. In comparison, the density functional 
values of Ref. |ll| of the spinodal pressure vary smoothly 
from the zero temperature value to the classical gas be- 
havior, having no flat region at low temperatures. The 
extrapolated experimental spinodal pressures of Hall and 
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cles form paths that span the width of the periodic cellc 
For a finite N, the superfiuid fraction should scale as, 
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FIG. 5. The temperature dependence of the spinodal 
density (•) and transition density (o) are shown with the 
experimental coexistence curve (dashed line) and superfiuid 
transition curve (solid line). The spinodal density exhibits 
a cusp-like feature as it intersects the superfiuid transition 
curve. The transition temperature for the hard sphere Bose 
gas of Ref. (A) and the ideal Bose gas (long dashed line) 
are shown for comparison. 



Marisli3 were limited to temperatures above 2.2 K or at 
K. 

It is illustrative to plot the liquid spinodal density ver- 
sus temperature, as is shown in Fig. |5|. The liquid spin- 
odal density shows a temperature dependence similar to 
the experimental liquid coexistence density (dashed line) : 
a temperature independent region at low temperatures 
with a slight cusp like behavior near where the experi- 
mental transition line (solid line) approaches. The transi- 
tion line at negative pressures as determined from PIMC 
simulation is discussed below. 



B. Superfiuid Transition 



Ps(T) 



L- x Q{L x l v t) 



(6) 



where t — (T — T\)/T\ is the reduced temperature. For 
the critical exponent of the correlation length, v, experi- 
ment gives v = 0.67. By assuming that the scaling func- 
tion Q is linear with respect to its argument near T = T\, 
the parameters v and T\ can be varied to minimize the 
distance of the data to Q. The universal constant Q{CfU 
see Table ||, has been calculated for the 3D XY modeled 
as 0.49(1). 

Using systems of 16, 32 and 64 particles at a fixed 
density, the superfiuid fraction was determined at sev- 
eral temperatures and scaled according to Eq. ^. Once 
the transition temperature was found, the corresponding 
pressure was found by interpolating the TV = 64 pressure 
data. The superfiuid transition values plotted in Fig. || 
and Fig. || are listed in Table 0. The superfiuid transition 
data fall along the experimental superfiuid transition line 
if it extended linearly to negative pressures. It appears 
that the transition line intersects the spinodal line at a 
temperature of 2.2 K. 

For comparispa, the transition temperature for a hard 
sphere Bose gasLJ (a system which does not have a liquid- 
gas transition) and an ideal Bose gas are shown in Fig. ||. 
As was explained by GrutertJ in a study of superfluidity 
for a hard sphere Bose gas , at low densities spatial fluc- 
tuations are important and clusters are likely to form, 
which inhibit macroscopic exchange cycles. At moder- 
ate densities the system is more homogeneous, allowing 
the exchange cycles necessary for superfluidity to form 
and the transition temperature becomes greater than the 
ideal Bose gas value by 7 %. At high densities, exchange 
is once again inhibited due to an increase in effective mass 
of the particles lowering the transition temperature be- 
low the ideal Bose gas value. Note that in comparing the 
temperature dependence of the transition densities, the 
pressures of the ideal and hard sphere Bose gases are pos- 
itive while the pressure of the 4 He transition goes from 
positive to negative as the density decreases. 



It is not clear what is the temperature dependencei-of 
the superfiuid transition line at negative pressuresl!3ll3. 
Experimentally, at positive pressures, as the density is 
reduced, the transition temperature increases. Further 
reducing the density, the transition line could extend 
into the negative pressure region. Additionally, along 
an isotherm, as the density is reduced to the spinodal 
density, .superfluidity will vanish as the speed of sound 
vanishesliZI. In order to accurately determine the super- 
fluid transition temperature we used finite size scaling of 
the superfiuid fractior£3. The direct estimator for the 
superfiuid fraction uses the winding number, which is a 
measure of the degree to which the particle exchange cy- 



C. Isothermal Sound Velocity 

The scaling of experimental isothermal sound veloc- 
ity data to pressure is controversiaL-, Early thermody- 
namic arguments put forth by Marisc^l indicated that at 



TABLE II. Superfiuid transition values. 



Tx (K) 


Pa (g/cc) 


P\ (bars) 


V 


Q(0) a 


2.20(2) 


0.11463 


-7.014(9) 


0.63(4) 


0.49(2) 


2.21(2) 


0.12669 


-6.034(11) 


0.65(4) 


0.55(2) 



Q(0) = h 2 pl /3 Q(0)/m 4 k B T x 
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FIG. 6. The isothermal speed of sound from Eq. is 
plottecLversus pressure for varying temperature. It has been 
shownEj that c 3 should vary linearly with pressure, vanishing 
at the spinodal pressure. The straight lines are fits to the 
PIMC data (filled symbols) and are to guide the eye. The 
experimental data collected in Ref. [l^ are indicated by the 
open symbols and are limited to positive pressures. 
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FIG. 7. The excitation spectra from PIMC simulations 
at a temperature of 1.0 K and densities of 0.114 ffAx and 
0.145 g/cc are compared to the experimental curveEa (solid 
line). At a density of 0.114 g/cc, the system is nearly at 
the spinodal pressure. Tangents at k — are shown (dashed 
lines), where the velocity is determined by from Eq. M. The 
lowest value of k for the simulations is determined by the size 
of the computational cell. 



T = K the speed of sound should scale with pressure 
as c 4 ~ p — p sl where .P.; is the liquid spinodal pressure. 
Maris' later analysisEZlEj indicated that this was not cor- 
rect in the pressure range accessible to experiment. A 
renormalization group analysis of the problem by MarisEZI 
yielded a scaling relation of c 3 ~ P — P s i for pressures 
near P s i, which proved to be very representative of the 
data of both 4 He and 3 Bej r -,There is some disagreement 
with this interpretatiorjcjfij. Accordingly, the c 4 scal- 
ing should occur only within a short distance from P s i 
while the c 3 scaling would be recovered at experimen- 
tally observed values of pressure. However, this seems 
counter to the renormalization approach where it would 
be expected that the exponent should be 4 away from the 
critical pressure and 3 some small distance from P s i. At 
issue is whether or not the second derivative of pressure 
with respect to density is identically zero at the spinodal 
density. Recent microscopic calculations by Campbell 
et al.13 show a different power law when very near the 
spinodal point. While the above discussionJs limited to 
T = K, the c 3 scaling behavior is observedE3 for all tem- 
peratures up to the critical temperature of T = 5.2 K. 



In Fig. ^ the isothermal sound velocity from PIMC 
simulation is plotted as a function of pressure, being de- 



termined numerically from from the PIMC pressure data, 

(7) 



TO4C 



dP 
dp ' 



As can be seen in the figure, the PIMC data is in agree- 
ment with the relation c 3 ~ P — P s i although there is 
some numerical error in taking the derivative. The ex- 
perimental data exhibits a smaller slope than the simula- 
tion data but the spinodal pressures (ordinate intercept) 
are in agreement. 



D. Excitation Spectrum 

The isothermal speed of sound is, seen to vanish on the 
spinodal line. According to MarisEj, in the quasi-particle 
picture of liquid helium, a transformation of the phonon 
branch to a free particle-like dispersion is expected to oc- 
cur as the liquid spinodal line is neared. Additionally, at 
negative pressures, a reduction of the maxon peak may 
occur in the excitation spectrum. The effect of lower- 
ing the maxon peak is a spinodal pressure which has a 
temperature dependence similar to that in Fig. || 

At a temperature of 1.0 K the excitation spectrum 
at two densitiesJaas been calculated using a maximum 
entropy analysisEl on the PIMC intermediate scattering 
function data and is shown in Fig. ^. Estimates for the 
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uncertainties of these values are difficult to make. The 
dynamic structure function is the Fourier transform of 
the intermediate scattering function, which is analyt- 
ically continued to imaginary time. In this instance, 
the transform is numerically underdctcrmincd with no 
unique solution. The peaks in the calculated dynamic 
structure function reproduce the experimental spectrum 
but the calculated widths of the dynamic-structure func- 
tion are too large in the superfluid phaseEI The neutron 
scattering data at T j^, 0.35 K and density 0.145 g/cc, 
compiled by DonnellycJ, is shown for comparison. The 
PIMC data at a density of 0.145 g/cc and pressure 
~ —0.603 bar agrees with the neutron data. At a den- 
sity of 0.112 g/cc and pressure ~ —7.72 bar (very nearly 
the liquid spinodal density) , there is an obvious decrease 
in the maxon peak but since the maxon peak is density 
dependent, it is not clear whether the effect is solely the 
result of the proximity to the spinodal. The phonon part 
of the spectrum is at momenta less then the lowest wave 
vector present in the system being simulated. However, 
the slope of the spectrum at zero wavevector is known, 
Eq. |7[ and is shown (dashed lines) for each value of the 
density. 



IV. CONCLUSIONS 

The spinodal line, the superfluid transition and the ex- 
citation spectrum for liquid helium at negative pressures 
have been investigated using PIMC. The spinodal line at 
low temperatures shows substantial agreement with the 
density functional calculations of Ref. [tl] as well as at 
temperatures above the superfluid transition. The calcu- 
lated temperature independence at low temperatures and 
the sudden increase in pressure as the superfluid transi- 
tion is neared is consistent with the conjectures of Ref. [l^. 
At higher temperatures there is general good agreement 
of the spinodal pressure with all calculations in the nor- 
mal fluid phase. The superfluid transition is found to be 
consistent with the experimental transition values and 
approaches the spinodal line near 2.2 K. Finally, the ex- 
citation spectrum exhibits changes in the maxon area 
consistent with the low temperature behavior of the liq- 
uid spinodal line. 
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